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Abstract 

It is known that scalar-tensor theory of gravity admits regular crossing of the phantom divide 
line wde = — 1 for dark energy, and existing viable models of present dark energy for its particular 
case ~ fiR) gravity - possess one such crossing in the recent past, after the end of the matter 
dominated stage. It was recently noted that during the future evolution of these models the dark 
energy equation of state wde may oscillate with an arbitrary number of phantom divide crossings. 
In this paper we prove that the number of crossings can be infinite, present an analytical condition 
for the existence of this effect and investigate it numerically. With the increase of the present 
mass of the scalaron (a scalar particle appearing in f{R) gravity) beyond the boundary of the 
appearance of such oscillations, their amplitude is shown to decrease very fast. As a result, the 
effect quickly becomes small and its beginning is shifted to the remote future. 
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I. INTRODUCTION 



The accelerating expansion of the present Universe is confirmed by current precise ob- 
servational data such as type la supernovae [1, 2], anisotropy of cosmic microwave back- 
ground [3], large scale structure [4] and baryon acoustic oscillations [5, 6]. The standard 
cosmological constant (A)-Cold-Dark-Matter (CDM) model is indeed able to explain these 
observational results within observational errors. In this model a cosmological constant is 
regarded as a new fundamental physical constant. However, the required value of the cos- 
mological constant is very tiny compared with any known physical scales. Thus, its relation 
to the standard quantum theory of known particles and fields is not understood today al- 
though some nonperturbative effects may naturally generate such a small quantity, see e.g. 
Refs. [7, 8] . More generally, a source of the current cosmic acceleration is called dark energy 
(DE). Further we shall use the more detailed term "present DE" for it to distinguish it from 
primordial DE which was responsible for another accelerated expansion regime, dubbed in- 
flation, which occurred in the very early Universe [9-11]. The relation of primordial DE to 
the known elementary particles has not been established, too. 

Both primordial and present DE can have either a physical origin (some new physical 
fields of matter) or a geometrical one. In the latter case, the Einstein gravity becomes 
modified. One of the simplest and self-consistent generalizations of the Einstein gravity is 
/(i?) gravity which incorporates a new phenomcnological function f{R) of the Ricci scalar 
R (with (ff/dR^ not identically zero) into the action, see Eq. (3) below. For a long time this 
theory of gravity was known to contain viable infiationary models, among them the simplest 
one introduced in Ref. [9] that remains in agreement with the most recent observational 
data. Thus, f{R) gravity can successively describe primordial DE. Rather recently, after 
many unsuccessful trials, viable models of present dark energy were found [12-14] which 
provide non-trivial alternatives to the standard ACDM model. ^ This theory is a special 
class of the scalar-tensor theory of gravity with the vanishing Brans-Dicke parameter uj Bo- 
lt contains a new scalar degree of freedom dubbed "scalaron" in Ref. [9], thus, it is a 

^ In order not to destroy the standard early Universe cosmology, including the recombination, the correct 
Big Bang nucleosynthesis and inflation of any kind, these models of present DE possessing a non-trivial 
form of f{R) in the low-i?, i? > region have to be further generalized by changing the behaviour of f{R) 
at large R and by extending it to the region of negative R, see Ref. [15] . However, this generalization is 
not important for our study. 
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nonperturbative generalization of the Einstein gravity. Prom the quantum point of view, 
scalaron is a massive spin-0 particle which mass depends on R. We consider f{R) gravity as 
a phenomenological macroscopic theory of gravity, alternative to the Einstein one, without 
discussing its microscopic origin.^ 

The existence of this additional degree of freedom imposes a number of constraints on 
the functional form of f{R) in viable cosmological models. In particular, in order to have 
the correct Newtonian limit, as well as the standard matter-dominated stage with the scale 
factor behaviour a{t) oc t^^^ driven by cold dark matter and baryons, the following conditions 
should be fulfilled for R^ Rq where Rq = R{to) ^ Hq, to is the present moment and Hq is 
the Hubble constant, and up to curvatures in the centre of neutron stars: 

\fiR)-R\<.R, Rf"{R)<^l. (1) 

Here the prime denotes the derivative with respect to the argument R. Furthermore, f{R) 
should satisfy the following conditions to guarantee both that Newtonian gravity solutions 
are stable and that the standard matter-dominated Friedmann-Robertson- Walker (FRW) 
stage remains an attractor with respect to an open set of neighbouring generic cosmological 
solutions in f{R) gravity: 

f{R) > 0, fiR) > 0. (2) 

In quantum language, the first condition means that gravity is attractive and graviton is not 
a ghost, while the second one - that scalaron is not a tachyon. Specific functional forms of 
f{R) satisfying these conditions, as well as laboratory and Solar system tests of gravity, and 
possessing a future stable (or at least metastable) de Sitter stage that is required for correct 
description of observable properties of present DE, have been proposed in Refs. [12-14], and 
much work has been carried out on their cosmological consequences. 

In order to describe the difference between FRW background solutions of f{R) gravity 
and the ACDM model, it is useful to introduce the effective equation-of-state (EoS) pa- 
rameter for DE wde = -Pde/pde where the effective pressure Pde and the effective energy 

^ Note the simplest possible mechanism that has attracted a new interest recently: a scalar field (/) with 
some potential and the non- minimal coupling — ^i?</)^/2 to the Einstein gravity in the limit of a very large 
negative ^ (i.e. the sign of coupling is opposite to that of the conformally coupled case). However, this 
mechanism leads to df/dR > 1. So, while sufiicient to produce the functional form of f{R) needed for 
successful inflationary models, it is not useful for construction of viable models of present DE. 
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density pde of DE are determined using the Einsteinian representation of gravitational field 
equations, see Eqs. (9), (10) below. Another independent parameter which describes scalar 
(density) perturbations on a FRW background is the gravitational growth index 7 defined 
as d\n5 /d\na = Qmiz)"'''^^ where 5 = 5pm/ Pm and Vt^ = STrGprn/SH"^ are a matter density 
fluctuation and the density parameter for matter, respectively. In f{R) gravity, Wde is time 
dependent and 7 is time and scale dependent whilst they keep the constant value wde = — 1 
and 7 6/11 in the ACDM model. Time and scale dependency of 7 generate an addi- 
tional transfer function for matter density fluctuation that constrains the model parameter 
region [16-18]. 

One of the most interesting features of geometrical DE distinguishing it from physical 
DE based on non-ghost physical fields minimally coupled to gravity, like quintessence, is 
the possibility of phantom behaviour, wbe < —1, of DE. Moreover, this behaviour may 
well be temporary with DE becoming normal, wde > —1, after smooth crossing of the 
phantom boundary wbe = — 1- In particular, models of geometrical DE based on scalar- 
tensor gravity were long known to admit this property [19]. f{R) gravity is a particular 
case of scalar-tensor gravity, so it permits phantom behaviour of DE and smooth crossing 
of the phantom boundary, too. Existing observational data do not exclude the possibility of 
phantom behaviour of DE (though they do not specifically favour it, too) for the following 
simple reason: as was noted above, DE in the particular form of an exact cosmological 
constant A is in a good agreement with all data. But since w\ = —1, it lies exactly at the 
phantom boundary. Thus, any small deviation of DE from A to the direction of decreasing 
wde results in its phantom behaviour. So, theory has to be prepared for this possibility 
that explains large interest in DE models admitting it. However, it is clear already that 
this "phantomness" should be small. In particular, if it is assumed for simplicity that 
wbe — const, when |wde + 1| < 0.1 at the approximately 2a confidence level [3]. 

Moreover, viable f{R) models of present DE [12-14] generically exhibit phantom be- 
haviour during the matter-dominated stage and one recent crossing of the phantom divide 
Wde = —1 even in the case of the smoothest behaviour of a FRW scale factor a{t), when 
there were no superimposed small oscillations of a(t) in the past (in quantum language, no 
condensate of primordial scalarons with the zero momentum) [12, 15, 16, 20]. From the 
physical point of view, the absence of primordial scalarons in the viable f{R) models of 
present DE is needed in order not to destroy the standard early Universe cosmology and 
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it can be achieved by primordial inflation of any kind, see Ref. [15] for a detailed discus- 
sion. Using numerical calculations, it has been recently shown that even in this smoothest 
case the EoS parameter wde can oscillate around the future de Sitter solution in these DE 
models [21], see also Ref. [22]. 

To investigate the phenomenon of multiple crossing of the phantom divide in more detail 
and analytically, in the present paper we prove that this crossing can indeed occur infinitely 
many times during the future evolution of viable f{R) models of present DE if the scalaron 
mass at a future stable de Sitter stage in these models is sufficiently large. Though this 
phenomenon is not directly observable since it refers to remote future, it is interesting and 
important from the theoretical point of view. Also, it is possible to check from observational 
data at the present moment if the derived analytical criterion for the existence of the infinite 
number of crossings is satisfied or not. 

Thus, the present paper focuses on the oscillatory behaviour of wde around the phantom 
divide wde = —1 in the future. In Sec. II, we review the stability conditions and the 
condition of the existence of a stable future de Sitter stage in f{R) gravity, and derive the 
condition for the existence of the infinite number of these oscillations analytically using 
the perturbation theory around the dc Sitter solution. In Sec. Ill, we focus on the specific 
viable model of present DE in f{R) gravity and present results from numerical calculations 
relating this condition to observable properties of the Universe at the present time. Sec. IV 
is devoted to conclusions and discussion. 

II. THE CRITERIA 

The action studied is of the form 



where f{R) is a function of the Ricci scalar R and Sm denotes matter action with the 
minimal coupling to gravity. Field equations are derived as 




(3) 



^T^GT^^ = (1 - F)R,, f)9,u + (VmV. - g^.n)F 



(4) 



(5) 
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where F — df/dR. We use the representation (4), (5) to define the effective energy- 
momentum tensor of DE. The (0, 0) and (i, i) components of the field equations are 

BF — f 

3FH^ = — ^ - 3HF + SttGp, (6) 

6F- = RF - f -3(F + HF) -8nG(p + 3P). (7) 
a 

It is also useful to use the trace equation: 

RF-2f + SnF = SttGT. (8) 

Thus, the effective energy density, the pressure and the EoS parameter of DE have the form: 

SnGpoE = 3//^ - SnGp = -3(1 - F) - + ^—^ - 3HF, (9) 

CI 2 

SttCPde = -^H - - SnGP = (1 - F) ( - + 2hA - ^—^ + F + 2HF, (10) 



2{1- F){-a/a + H') + F-HF 
wde + 1 = -■ (11) 

-3{1- Fyd/a + {R- f )/2-3HF ^ ' 

During an asymptotic de Sitter regime, the matter density decreases rapidly as p oc e~^^^* 
and soon can be neglected. Therefore, it follows from Eq. (8) that a constant value of the 
Ricci scalar R — R^ — const = 12H^ at a de Sitter regime is given by a root of the algebraic 
equation 

2/i = RiF, (12) 

where /i = f{Ri) and Fi = F{Ri). At the de Sitter regime, DE is characterized by 

SttGpde,! = -SttCPde,! = X' ^^"^ ^DE,i = -1- 

To investigate the stabihty of the future de Sitter solution and the possibihty of oscillatory 
behaviour around it, we expand Eqs. (6), (8) in the perturbation series with respect to 
SR = R — Ri and 5H = H — Hi. In the first order in 5R and 5H, 

where prime denotes the derivative with respect to the number of e- folds A'" = Ina = 
— ln(l + 2;) and Fri = Fr(-Ri) = dF{Ri)/dR. Although the matter density term in the 
right-hand side has the zero order, we include it because pm oc e~^^^* is much smaller than 
background quantities at the future de Sitter stage. 
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SR" + 3SR' +-±-^{^-R,)SR= ^^^P^ , (14) 



Eq. (14) is solved as a sum of the homogeneous solution of Eq. (14) with the zero 
right-hand side, SRosc, and the special solution of the non-homogeneous equation SR^ec' 

5R = SRdec + SRosc- (15) 
Since Pm = Pmo^"^^, SR^ec is obtained as 

SR.^ = '"'t: e^'"- (16) 

Pi — KiPri 

Notice that it describes a monotonically decaying mode. 

On the other hand, the homogeneous solution 5i?osc may have decaying, growing and 
oscillatory behaviour. In the case of monotonic behaviour (both roots of the characteristic 
equation for Eq. (14) are real), to keep the future dc Sitter solution stable (a stable node), 
the coefficient of the third term in the left-hand side of Eq. (14) should be positive. So, the 
stability condition is [23]: 

>Ri. (17) 



Fri 

The oscillatory behaviour occurs when the dc Sitter asymptote is a focus (complex roots). 
For this, the discriminant of the characteristic equation should be negative: 

> '^R. . (18) 



Fri 16 

Since the coefficient of the second term in the left-hand side of Eq. (14) is positive, the 

focus is always stable and the inequality (18) is stronger than (17). The criterion (18) of 
the oscillatory approach to the future de Sitter asymptote is equivalent to the condition 

' 3F^i 4 16 ' ^ ^ 

where Hi,Ri and Mi are the Hubble parameter, the scalar curvature and the scalaron mass 
at the future de Sitter state. If the oscillation condition is satisfied, 

SRosc = Ae-^^/^ sm{ujN + 0) (20) 

where ou = '^'\J r^Fri ~ it' ^ ^ integration constants. 

The perturbation of the EoS parameter 5wde = (<^-Pde + (^Pde)/pde,i is calculated from 
87rG'(pDE + Pde) = -2iir - SttGp^ and Eq. (13), 
4 



Swxi^ = — 
til 



3Fi 3 V Fi } \Fi 



(21) 



We decompose 5wde = Swdec + Swosc as 
4/1 



dw, 



dec 



Ri \Fi — RiFj 



Rl 



Ki 



1 j SttGp^oII + zf (22) 
cos{uN + <P) + 1 - l) sin{uN + 0)1 . (23) 



Thus, in the latter case of a stable de Sitter solution with oscillations, wde crosses the 
phantom boundary w^e — infinitely many times during the future evolution of the 
Universe. 

III. THE SPECIFIC MODEL 



We consider the following viable f{R) model [14]: 

f{R)=R + XR, + 



(24) 



where n and A are model parameters, and Rs is determined by the present observational 
data, namely, the ratio Rs/H^ is well fit by a simple power-law Rs/Hq = Cn\~^" with 
{n,Cn,Pn) = (2,4.16,0.953), (3,4.12,0.837), and (4,4.74,0.702), respectively [16]. 
From Eq. (12), the equation for de Sitter solutions is 

"1 + (n + l)r^ 



a{r) = r + 2A 



= 0, 



(25) 



_ (l+r2)"+i 

where r = Ri/Rg. It is obvious that the Minkowski space-time, r = 0, is one of the solutions. 
We denote the other positive solutions for Q;(r) = as = Ria/Rs and = Ru/Rg- 
Ta and Tft can be estimated by considering the limits of small and large r. For r <^ 1, 
a{r) ~ r[l — \n{n + l)r\ and for r » 1, a{r) ~ r — 2A. Therefore, the de Sitter solutions 
arc r = Ta — [\n{n + 1)]^"*^^^ and r = r^, c:^ 2A. Strictly speaking, these approximations are 
valid either for A ^ 1 or, in the case of Ta, for n ^ 1 (while A may be of the order of unity). 
However, it follows from the numerical analysis that the solutions for n = 2 and A = 3 are 
already close enough to these analytical estimations. 

One can check their stability and oscillatory behaviour using the stability parameter ^(r) 
and the oscillation parameter 7(r) which are derived from Eqs. (17) and (18): 

;i + r2)[(l + r2)"+i -2nAr] 



/3(r) 



2nA[(2n + l)r2 



r > 0, 



(l + r^)[(l + r^)"+^ -2nAr] _ 25 

2nA[(2n + l)r2- 1] 16^^ ' 



(26) 
(27) 
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TABLE I: Stable de Sitter solutions for various model parameters n and A. rf, = Rih/Rg is a stable 
de Sitter solution in terms of the normalized Ricci scalar. (3 and 7 denote the stability (26) and 
oscillation (27) parameters, respectively. 



n 


A 




Pin) 








2 


1 


1.64 


1.58 




6.56 X 


10-1 


2 


3 


5.99 


8.54 X 


10^ 


8.51 X 


102 


2 


10 


20.0 


3.23 X 


10^ 


3.23 X 


10^ 


3 


1 


1.94 


1.37 X 


10 


1.26 X 


10 


3 


3 


6.00 


1.53 X 


10^ 


1.53 X 


10^ 


3 


10 


20.0 


6.17 X 


10^ 


6.17 X 


10^ 


4 


1 


1.99 


5.07 X 


10 


4.95 X 


10 


4 


3 


6.00 


3.31 X 


10^ 


3.31 X 


10^ 


4 


10 


20.0 


1.44 X 


lO^o 


1.44 X 


lO^o 



Since 7(r) = /3(r) — 9r/16, there is no oscillations for an unstable de Sitter solution, as it 
should be. Prom these criteria we see that r — Va and r — are an unstable de Sitter 
solution and a stable de Sitter solution, respectively. The specific values are presented in 
the Table I. 

For a fixed n and various values of A, we obtain A^ and A^ as roots of /^(rt) = and 
j{i^b) = 0, respectively. Models are classified according to A being in the intervals A < 
XjS, Xjs < X < A-y, and A > A-y, and in each region a de Sitter solution r = r;, is unstable, 
stable without oscillations, and stable with oscillations, respectively. Although most of the 
parameters realize a stable de Sitter solution with oscillations (a stable focus), there exists 
a parameter region corresponding to a stable de Sitter solution without oscillation (a stable 
node). Fig. 1 suggests that such parameter regions are 0.944 < A < 0.970, 0.726 < A < 0.744 
and 0.608 < A < 0.622 for n = 2, 3 and 4, respectively. 

We integrate the evolution equations numerically. Initial condition are set at 2; = 10 using 
the ACDM model, and the present moment is determined by the condition = 0.27. Fig. 2 
shows that R approaches a stable de Sitter solution. It is seen from the right panel of Fig. 2 
that the perturbation theory with respect to 5R = R — Rn, is valid when z < —0.8 for 
n = 2, A = 1, and when z < —0.5 for n = 2, A = 3 or 10. The oscillation of 6R for 
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X 

FIG. 1: Values of the stability parameter /3 and the oscillation parameter 7 for stable de Sitter 
solutions for various model parameters. The parameter regions 7(rb) < < /3{rb) and 7(^5) > 
correspond to stable de Sitter solutions without oscillations and with oscillations, respectively. 

n = 2, A = 1 is clearly seen. For n = 2 and A = 3 or 10, oscillations exist, too, but their 
amplitude is so small that we cannot see them. To make them visible, we have subtracted 
the decaying mode SR^ec in Fig. 3. The analytic solution for 6Rosc fits the result well. 

Fig. 4 depicts the evolution of the EoS parameter for n = 2 and A = 1, 3, 10. The first 
phantom crossing occurred in the past at z ~ 0.5 in agreement with Ref. [16]. We subtract 
the decaying mode 6wdec and present the oscillation mode 6wosc in Fig. 5. The numerical 
results are very close to the analytic solutions for n = 2, A = 1 and 3. For n = 2, A = 10, 
the amplitude of the oscillations is small and the frequency is large, so that we cannot 
distinguish them from numerical noise. Finally, we present the case n = 2, A = 0.95 in 
Fig. 6 as an example of a non-oscillatory approach to the stable de Sitter solution. Note 
that the trajectories of 6R and 6w are convex upward and there is no oscillations indeed. 



IV. CONCLUSIONS AND DISCUSSION 



We have investigated conditions under which the effective EoS parameter wde of present 
DE in f{R) gravity can oscillate an infinite number of times around the phantom boundary 
wde = — 1 during the future evolution of the Universe. The analytical condition of the 
existence of this phenomenon, Eq. (18), is derived that depends on the properties of f{R) 
near a future stable de Sitter stage only. The physical sense of this condition is that the 
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FIG. 2: Future evolution of the Ricci scalar. It approaches the stable de Sitter solution which is 
presented in the Table I. 
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FIG. 3: Numerical results for 5R — SR^ec compared with the analytic solution for SRqsc- 



rest mass of the scalaron (a massive scalar particle which arises in f{R) gravity in addition 
to massless spin-2 graviton) should be sufficiently large at the future de Sitter stage. Thus, 
this phenomenon is generic. However, the amplitude of these oscillations has been shown to 
decrease fast with the increase of the scalaron mass beyond the boundary of the appearance 
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FIG. 4: Future evolution of the effective EoS parameter for dark energy. 




FIG. 5: Numerical results for (1 + wde) — <^u^dec compared with the analytic solution for 5wosc- 

of such oscillations. As a result, the effect quickly becomes small and its beginning is shifted 
to the remote future. For real scalaron masses lying below this boundary, the future stable 
de Sitter stage is reached without the phantom boundary crossing. Analytic solutions for 
the behaviour of Wbe near the phantom boundary have been obtained in the first order of 
the small quantity |wde + 1|- Generically they have a monotonically decaying part 6wdec 
and a damped harmonic oscillatory part Swosc- For a specific viable f{R) model of present 
DE energy, numerical integration of FRW background evolution has been performed which 
future behaviour is in a good agreement with the analytic formulas. 

All calculations have been done for the smoothest initial conditions in the past corre- 
sponding to the absence of a primordial homogeneous oscillating scalaron component. So, 
even in this case, an oscillating scalaron component (the condensate of scalarons with the 
zero momentum) arises around the present moment when the scalaron mass is comparable 
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FIG. 6: Future evolution of the effective EoS parameter for dark energy for n = 2 and A = 0.95. 
There is no oscillations. 



to the Hubble constant Hq (in the units where h = c= 1), and it quickly becomes dominant 
over the non-relativistic matter component (cold dark matter and baryons) in the future. 
But its effective energy-momentum tensor in turn soon becomes negligible compared to an 
effective cosmological constant producing the future stable de Sitter stage. For less smooth 
initial conditions, more phantom boundary crossings may occur in the past. But these ini- 
tial conditions are hardly compatible with the standard cosmology of the early Universe 
confirmed by numerous observational data. We hope to return to this question elsewhere. 

Finally, though the very phenomenon of multiple (and even infinite) number of phantom 
boundary crossings in the future is not directly observable, it is very interesting and impor- 
tant from the theoretical point of view. Also, as follows from our numerical calculations of 
the full evolution from the remote past to the remote future, the scalaron mass at the future 
de Sitter stage is close to its present value. Therefore, in principle it is possible to check 
from observational data describing the present and the past of our Universe if the derived 
analytical criterion for the existence of an infinite number of oscillations in is satisfied 
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or not. 
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